require(estimatr)

#### data preparation ####
FtF.data <- read.csv("FtF_data.csv")
online.data <- read.csv("online_data.csv")

# dichotomize education
FtF.data$college <- (FtF.data$education > 7) * 1
online.data$college <- (online.data$education > 7) * 1

# recategorize region
FtF.data$Visayas <- (FtF.data$region > 8 & FtF.data$region < 12) * 1
FtF.data$Mindanao <- (FtF.data$region > 11) * 1
FtF.data$area <- factor(ifelse(FtF.data$Visayas == 1, "Visayas", 
                               ifelse(FtF.data$Mindanao == 1, 
                                      "Mindanao", "Luzon")), 
                        levels = c("Luzon", "Visayas", "Mindanao"))
online.data$Visayas <- (online.data$region %in% c(2, 7, 10, 12, 16)) * 1
online.data$Mindanao <-(online.data$region %in% c(1, 5, 9, 14, 15, 17)) * 1
online.data$area <- factor(ifelse(online.data$Visayas == 1, "Visayas", 
                                  ifelse(online.data$Mindanao == 1, 
                                         "Mindanao", "Luzon")), 
                           levels = c("Luzon", "Visayas", "Mindanao"))

# dichotomize socio-economic class
FtF.data$lower <- (FtF.data$class > 3) * 1

# dichotomize ethnicity
FtF.data$nonTag <- (FtF.data$ethnicity > 1) * 1

# reorder internet use
FtF.data$internet.re <- factor(FtF.data$internet, levels = 5:0)

# reorder perceived neighborhood satisfaction with Duterte 
FtF.data$neighborhood.re <- factor(FtF.data$neighborhood, levels = 4:1)
online.data$neighborhood.re <- factor(online.data$neighborhood, levels = 4:1)

## function for subgroup comparison of SDB
SDB.group.test <- function(y, x, data, treat, block, DQ) {
  temp.data.1 <- subset(data, data[, x] == 0)
  n.obs.1 <- nrow(temp.data.1)
  list.reg.1 <- lm_robust(formula(paste0(y, " ~ ", treat)), data = temp.data.1, 
                          fixed_effects = paste0("~ ", block))
  list.est.1 <- as.numeric(list.reg.1$coefficients)
  list.se.1 <- as.numeric(list.reg.1$std.error)
  DQ.est.1 <- mean(temp.data.1[, DQ], na.rm = TRUE)
  DQ.se.1 <- sqrt(DQ.est.1 * (1 - DQ.est.1) / sum(! is.na(temp.data.1[, DQ])))
  SDB.est.1 <- DQ.est.1 - list.est.1
  SDB.se.1 <- sqrt(list.se.1 ^ 2 + DQ.se.1 ^ 2)
  SDB.p.1 <- (1 - pnorm(abs(SDB.est.1 / SDB.se.1))) * 2
  temp.data.2 <- subset(data, data[, x] == 1)
  n.obs.2 <- nrow(temp.data.2)
  list.reg.2 <- lm_robust(formula(paste0(y, " ~ ", treat)), data = temp.data.2, 
                          fixed_effects = paste0("~ ", block))
  list.est.2 <- as.numeric(list.reg.2$coefficients)
  list.se.2 <- as.numeric(list.reg.2$std.error)
  DQ.est.2 <- mean(temp.data.2[, DQ], na.rm = TRUE)
  DQ.se.2 <- sqrt(DQ.est.2 * (1 - DQ.est.2) / sum(! is.na(temp.data.2[, DQ])))
  SDB.est.2 <- DQ.est.2 - list.est.2
  SDB.se.2 <- sqrt(list.se.2 ^ 2 + DQ.se.2 ^ 2)
  SDB.p.2 <- (1 - pnorm(abs(SDB.est.2 / SDB.se.2))) * 2
  SDB.diff.est <- SDB.est.2 - SDB.est.1
  SDB.diff.se <- sqrt(list.se.1 ^ 2 + DQ.se.1 ^ 2 + list.se.2 ^ 2 + DQ.se.2 ^ 2)
  SDB.diff.p <- (1 - pnorm(abs(SDB.diff.est / SDB.diff.se))) * 2
  out <- list(summary = c(SDB.est.1, SDB.est.2, SDB.diff.est, 
                          SDB.diff.est + qnorm(0.025) * SDB.diff.se, 
                          SDB.diff.est + qnorm(0.975) * SDB.diff.se, 
                          SDB.diff.p), 
              details = list(SDB.est.1 = SDB.est.1, SDB.se.1 = SDB.se.1, 
                             SDB.p.1 = SDB.p.1, n.obs.1 = n.obs.1, 
                             SDB.est.2 = SDB.est.2, SDB.se.2 = SDB.se.2, 
                             SDB.p.2 = SDB.p.2, n.obs.2 = n.obs.2, 
                             SDB.diff.est = SDB.diff.est, 
                             SDB.diff.se = SDB.diff.se, 
                             SDB.diff.p = SDB.diff.p))
  names(out$summary) <- c("SDB est (x = 0)", "SDB est (x = 1)", "SDB diff est", 
                          "SDB diff lower", "SDB diff upper", "SDB diff p")
  out
}

#### significance tests of the differences in SDB ####
## gender
gender.result.FtF <- SDB.group.test("list", "female", FtF.data, 
                                    "treatment", "block", "DQ")
round(gender.result.FtF$summary, 3)

gender.result.online <- SDB.group.test("list", "female", online.data, 
                                       "treatment", "block", "DQ")
round(gender.result.online$summary, 3)

## education
edu.result.FtF <- SDB.group.test("list", "college", FtF.data, 
                                 "treatment", "block", "DQ")
round(edu.result.FtF$summary, 3)

edu.result.online <- SDB.group.test("list", "college", online.data, 
                                    "treatment", "block", "DQ")
round(edu.result.online$summary, 3)

## region
# comparison between Luzon and Visayas
LV.result.FtF <- SDB.group.test("list", "Visayas", 
                                subset(FtF.data, Mindanao != 1),
                                "treatment", "block", "DQ")
round(LV.result.FtF$summary, 3)

LV.result.online <- SDB.group.test("list", "Visayas", 
                                   subset(online.data, Mindanao != 1), 
                                   "treatment", "block", "DQ")
round(LV.result.online$summary, 3)

# comparison between Luzon and Mindanao
LM.result.FtF <- SDB.group.test("list", "Mindanao", 
                                subset(FtF.data, Visayas != 1),
                                "treatment", "block", "DQ")
round(LM.result.FtF$summary, 3)

LM.result.online <- SDB.group.test("list", "Mindanao", 
                                   subset(online.data, Visayas != 1), 
                                   "treatment", "block", "DQ")
round(LM.result.online$summary, 3)

# comparison between Visayas and Mindanao
VM.result.FtF <- SDB.group.test("list", "Mindanao", 
                                subset(FtF.data, Visayas == 1 | Mindanao == 1),
                                "treatment", "block", "DQ")
round(VM.result.FtF$summary, 3)

VM.result.online <- SDB.group.test("list", "Mindanao", 
                                   subset(online.data, Visayas == 1 | Mindanao == 1), 
                                   "treatment", "block", "DQ")
round(VM.result.online$summary, 3)

## socio-economic class
class.result.FtF <- SDB.group.test("list", "lower", FtF.data, 
                                   "treatment", "block", "DQ")
round(class.result.FtF$summary, 3)

#### Figure 4 ####
cairo_pdf("Figure_4.pdf", width = 5, height = 8, pointsize = 9)
layout(matrix(1:2, 2, 1))
par(mar = c(4, 0, 3, 1), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-1.95, 1), ylim = c(0, 19), 
     main = "Face-to-face survey", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(-0.6, c(1:3, 5:10, 12:14, 16:18), 1.05, c(1:3, 5:10, 12:14, 16:18), 
         lty = 3, col = "gray")
abline(v = seq(-0.5, 1, 0.25), col = "gray")
points(gender.result.FtF$details$SDB.est.1, 18, pch = 4)
segments(gender.result.FtF$details$SDB.est.1 + 
           qnorm(0.025) * gender.result.FtF$details$SDB.se.1, 18, 
         gender.result.FtF$details$SDB.est.1 + 
           qnorm(0.975) * gender.result.FtF$details$SDB.se.1, 18)
points(gender.result.FtF$details$SDB.est.2, 17, pch = 4)
segments(gender.result.FtF$details$SDB.est.2 + 
           qnorm(0.025) * gender.result.FtF$details$SDB.se.2, 17, 
         gender.result.FtF$details$SDB.est.2 + 
           qnorm(0.975) * gender.result.FtF$details$SDB.se.2, 17)
points(gender.result.FtF$details$SDB.diff.est, 16, pch = 19)
segments(gender.result.FtF$details$SDB.diff.est + 
           qnorm(0.025) * gender.result.FtF$details$SDB.diff.se, 16, 
         gender.result.FtF$details$SDB.diff.est + 
           qnorm(0.975) * gender.result.FtF$details$SDB.diff.se, 16)
text(gender.result.FtF$details$SDB.diff.est, 16, 
     sprintf("%5.3f", gender.result.FtF$details$SDB.diff.est), 
     cex = 0.6, pos = 1)
points(edu.result.FtF$details$SDB.est.1, 14, pch = 4)
segments(edu.result.FtF$details$SDB.est.1 + 
           qnorm(0.025) * edu.result.FtF$details$SDB.se.1, 14, 
         edu.result.FtF$details$SDB.est.1 + 
           qnorm(0.975) * edu.result.FtF$details$SDB.se.1, 14)
points(edu.result.FtF$details$SDB.est.2, 13, pch = 4)
segments(edu.result.FtF$details$SDB.est.2 + 
           qnorm(0.025) * edu.result.FtF$details$SDB.se.2, 13, 
         edu.result.FtF$details$SDB.est.2 + 
           qnorm(0.975) * edu.result.FtF$details$SDB.se.2, 13)
points(edu.result.FtF$details$SDB.diff.est, 12, pch = 19)
segments(edu.result.FtF$details$SDB.diff.est + 
           qnorm(0.025) * edu.result.FtF$details$SDB.diff.se, 12, 
         edu.result.FtF$details$SDB.diff.est + 
           qnorm(0.975) * edu.result.FtF$details$SDB.diff.se, 12)
text(edu.result.FtF$details$SDB.diff.est, 12, 
     sprintf("%5.3f", edu.result.FtF$details$SDB.diff.est), 
     cex = 0.6, pos = 1)
points(LV.result.FtF$details$SDB.est.1, 10, pch = 4)
segments(LV.result.FtF$details$SDB.est.1 + 
           qnorm(0.025) * LV.result.FtF$details$SDB.se.1, 10, 
         LV.result.FtF$details$SDB.est.1 + 
           qnorm(0.975) * LV.result.FtF$details$SDB.se.1, 10)
points(LV.result.FtF$details$SDB.est.2, 9, pch = 4)
segments(LV.result.FtF$details$SDB.est.2 + 
           qnorm(0.025) * LV.result.FtF$details$SDB.se.2, 9, 
         LV.result.FtF$details$SDB.est.2 + 
           qnorm(0.975) * LV.result.FtF$details$SDB.se.2, 9)
points(VM.result.FtF$details$SDB.est.2, 8, pch = 4)
segments(VM.result.FtF$details$SDB.est.2 + 
           qnorm(0.025) * VM.result.FtF$details$SDB.se.2, 8, 
         VM.result.FtF$details$SDB.est.2 + 
           qnorm(0.975) * VM.result.FtF$details$SDB.se.2, 8)
points(LV.result.FtF$details$SDB.diff.est, 7, pch = 19)
segments(LV.result.FtF$details$SDB.diff.est + 
           qnorm(0.025) * LV.result.FtF$details$SDB.diff.se, 7, 
         LV.result.FtF$details$SDB.diff.est + 
           qnorm(0.975) * LV.result.FtF$details$SDB.diff.se, 7)
text(LV.result.FtF$details$SDB.diff.est, 7, 
     sprintf("%5.3f", LV.result.FtF$details$SDB.diff.est), 
     cex = 0.6, pos = 1)
points(LM.result.FtF$details$SDB.diff.est, 6, pch = 19)
segments(LM.result.FtF$details$SDB.diff.est + 
           qnorm(0.025) * LM.result.FtF$details$SDB.diff.se, 6, 
         LM.result.FtF$details$SDB.diff.est + 
           qnorm(0.975) * LM.result.FtF$details$SDB.diff.se, 6)
text(LM.result.FtF$details$SDB.diff.est, 6, 
     sprintf("%5.3f", LM.result.FtF$details$SDB.diff.est), 
     cex = 0.6, pos = 1)
points(VM.result.FtF$details$SDB.diff.est, 5, pch = 19)
segments(VM.result.FtF$details$SDB.diff.est + 
           qnorm(0.025) * VM.result.FtF$details$SDB.diff.se, 5, 
         VM.result.FtF$details$SDB.diff.est + 
           qnorm(0.975) * VM.result.FtF$details$SDB.diff.se, 5)
text(VM.result.FtF$details$SDB.diff.est, 5, 
     sprintf("%5.3f", VM.result.FtF$details$SDB.diff.est), 
     cex = 0.6, pos = 1)
points(class.result.FtF$details$SDB.est.1, 3, pch = 4)
segments(class.result.FtF$details$SDB.est.1 + 
           qnorm(0.025) * class.result.FtF$details$SDB.se.1, 3, 
         class.result.FtF$details$SDB.est.1 + 
           qnorm(0.975) * class.result.FtF$details$SDB.se.1, 3)
points(class.result.FtF$details$SDB.est.2, 2, pch = 4)
segments(class.result.FtF$details$SDB.est.2 + 
           qnorm(0.025) * class.result.FtF$details$SDB.se.2, 2, 
         class.result.FtF$details$SDB.est.2 + 
           qnorm(0.975) * class.result.FtF$details$SDB.se.2, 2)
points(class.result.FtF$details$SDB.diff.est, 1, pch = 19)
segments(class.result.FtF$details$SDB.diff.est + 
           qnorm(0.025) * class.result.FtF$details$SDB.diff.se, 1, 
         class.result.FtF$details$SDB.diff.est + 
           qnorm(0.975) * class.result.FtF$details$SDB.diff.se, 1)
text(class.result.FtF$details$SDB.diff.est, 1, 
     sprintf("%5.3f", class.result.FtF$details$SDB.diff.est), 
     cex = 0.6, pos = 1)
text(-1.95, 19, "Gender", pos = 4, font = 2)
text(-1.8, 18, paste0("Male (N=", 
                      format(gender.result.FtF$details$n.obs.1, big.mark = ","), 
                      ")"), pos = 4)
text(-1.8, 17, paste0("Female (N=", 
                      format(gender.result.FtF$details$n.obs.2, big.mark = ","), 
                      ")"), pos = 4)
text(-1.8, 16, "Difference b/w groups", pos = 4)
text(-1.95, 15, "Education", pos = 4, font = 2)
text(-1.8, 14, paste0("Lower than college (N=", 
                      format(edu.result.FtF$details$n.obs.1, big.mark = ","), 
                      ")"), pos = 4)
text(-1.8, 13, paste0("College or higher (N=", 
                      format(edu.result.FtF$details$n.obs.2, big.mark = ","), 
                      ")"), pos = 4)
text(-1.8, 12, "Difference b/w groups", pos = 4)
text(-1.95, 11, "Region", pos = 4, font = 2)
text(-1.8, 10, paste0("Luzon (N=", 
                      format(LV.result.FtF$details$n.obs.1, big.mark = ","), 
                      ")"), pos = 4)
text(-1.8, 9, paste0("Visayas (N=", 
                     format(LV.result.FtF$details$n.obs.2, big.mark = ","), 
                     ")"), pos = 4)
text(-1.8, 8, paste0("Mindanao (N=", 
                     format(VM.result.FtF$details$n.obs.2, big.mark = ","), 
                     ")"), pos = 4)
text(-1.8, 7, "Difference b/w L & V", pos = 4)
text(-1.8, 6, "Difference b/w L & M", pos = 4)
text(-1.8, 5, "Difference b/w V & M", pos = 4)
text(-1.95, 4, "Socio-economic class", pos = 4, font = 2)
text(-1.8, 3, paste0("Upper class (N=", 
                     format(class.result.FtF$details$n.obs.1, big.mark = ","), 
                     ")"), pos = 4)
text(-1.8, 2, paste0("Lower class (N=", 
                     format(class.result.FtF$details$n.obs.2, big.mark = ","), 
                     ")"), pos = 4)
text(-1.8, 1, "Difference b/w groups", pos = 4)
axis(1, seq(-0.5, 1, 0.25), lwd = 0.5, labels = NA)
mtext(sprintf("%4.2f", seq(-0.5, 1, 0.25)), 
      side = 1, line = 1, at = seq(-0.5, 1, 0.25))
mtext("Magnitude of SDB (or its difference)", side = 1, line = 2.5, at = 0.25)
plot(NULL, NULL, bty = "n", xlim = c(-1.95, 1), ylim = c(0, 19), 
     main = "Online survey", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(-0.6, c(1:3, 5:10, 12:14, 16:18), 1.05, c(1:3, 5:10, 12:14, 16:18), 
         lty = 3, col = "gray")
abline(v = seq(-0.5, 1, 0.25), col = "gray")
points(gender.result.online$details$SDB.est.1, 18, pch = 4)
segments(gender.result.online$details$SDB.est.1 + 
           qnorm(0.025) * gender.result.online$details$SDB.se.1, 18, 
         gender.result.online$details$SDB.est.1 + 
           qnorm(0.975) * gender.result.online$details$SDB.se.1, 18)
points(gender.result.online$details$SDB.est.2, 17, pch = 4)
segments(gender.result.online$details$SDB.est.2 + 
           qnorm(0.025) * gender.result.online$details$SDB.se.2, 17, 
         gender.result.online$details$SDB.est.2 + 
           qnorm(0.975) * gender.result.online$details$SDB.se.2, 17)
points(gender.result.online$details$SDB.diff.est, 16, pch = 19)
segments(gender.result.online$details$SDB.diff.est + 
           qnorm(0.025) * gender.result.online$details$SDB.diff.se, 16, 
         gender.result.online$details$SDB.diff.est + 
           qnorm(0.975) * gender.result.online$details$SDB.diff.se, 16)
text(gender.result.online$details$SDB.diff.est, 16, 
     sprintf("%5.3f", gender.result.online$details$SDB.diff.est), 
     cex = 0.6, pos = 1)
points(edu.result.online$details$SDB.est.1, 14, pch = 4)
segments(edu.result.online$details$SDB.est.1 + 
           qnorm(0.025) * edu.result.online$details$SDB.se.1, 14, 
         edu.result.online$details$SDB.est.1 + 
           qnorm(0.975) * edu.result.online$details$SDB.se.1, 14)
points(edu.result.online$details$SDB.est.2, 13, pch = 4)
segments(edu.result.online$details$SDB.est.2 + 
           qnorm(0.025) * edu.result.online$details$SDB.se.2, 13, 
         edu.result.online$details$SDB.est.2 + 
           qnorm(0.975) * edu.result.online$details$SDB.se.2, 13)
points(edu.result.online$details$SDB.diff.est, 12, pch = 19)
segments(edu.result.online$details$SDB.diff.est + 
           qnorm(0.025) * edu.result.online$details$SDB.diff.se, 12, 
         edu.result.online$details$SDB.diff.est + 
           qnorm(0.975) * edu.result.online$details$SDB.diff.se, 12)
text(edu.result.online$details$SDB.diff.est, 12, 
     sprintf("%5.3f", edu.result.online$details$SDB.diff.est), 
     cex = 0.6, pos = 1)
points(LV.result.online$details$SDB.est.1, 10, pch = 4)
segments(LV.result.online$details$SDB.est.1 + 
           qnorm(0.025) * LV.result.online$details$SDB.se.1, 10, 
         LV.result.online$details$SDB.est.1 + 
           qnorm(0.975) * LV.result.online$details$SDB.se.1, 10)
points(LV.result.online$details$SDB.est.2, 9, pch = 4)
segments(LV.result.online$details$SDB.est.2 + 
           qnorm(0.025) * LV.result.online$details$SDB.se.2, 9, 
         LV.result.online$details$SDB.est.2 + 
           qnorm(0.975) * LV.result.online$details$SDB.se.2, 9)
points(VM.result.online$details$SDB.est.2, 8, pch = 4)
segments(VM.result.online$details$SDB.est.2 + 
           qnorm(0.025) * VM.result.online$details$SDB.se.2, 8, 
         VM.result.online$details$SDB.est.2 + 
           qnorm(0.975) * VM.result.online$details$SDB.se.2, 8)
points(LV.result.online$details$SDB.diff.est, 7, pch = 19)
segments(LV.result.online$details$SDB.diff.est + 
           qnorm(0.025) * LV.result.online$details$SDB.diff.se, 7, 
         LV.result.online$details$SDB.diff.est + 
           qnorm(0.975) * LV.result.online$details$SDB.diff.se, 7)
text(LV.result.online$details$SDB.diff.est, 7, 
     sprintf("%5.3f", LV.result.online$details$SDB.diff.est), 
     cex = 0.6, pos = 1)
points(LM.result.online$details$SDB.diff.est, 6, pch = 19)
segments(LM.result.online$details$SDB.diff.est + 
           qnorm(0.025) * LM.result.online$details$SDB.diff.se, 6, 
         LM.result.online$details$SDB.diff.est + 
           qnorm(0.975) * LM.result.online$details$SDB.diff.se, 6)
text(LM.result.online$details$SDB.diff.est, 6, 
     sprintf("%5.3f", LM.result.online$details$SDB.diff.est), 
     cex = 0.6, pos = 1)
points(VM.result.online$details$SDB.diff.est, 5, pch = 19)
segments(VM.result.online$details$SDB.diff.est + 
           qnorm(0.025) * VM.result.online$details$SDB.diff.se, 5, 
         VM.result.online$details$SDB.diff.est + 
           qnorm(0.975) * VM.result.online$details$SDB.diff.se, 5)
text(VM.result.online$details$SDB.diff.est, 5, 
     sprintf("%5.3f", VM.result.online$details$SDB.diff.est), 
     cex = 0.6, pos = 1)
text(-1.95, 19, "Gender", pos = 4, font = 2)
text(-1.8, 18, paste0("Male (N=", 
                      format(gender.result.online$details$n.obs.1, big.mark = ","), 
                      ")"), pos = 4)
text(-1.8, 17, paste0("Female (N=", 
                      format(gender.result.online$details$n.obs.2, big.mark = ","), 
                      ")"), pos = 4)
text(-1.8, 16, "Difference b/w groups", pos = 4)
text(-1.95, 15, "Education", pos = 4, font = 2)
text(-1.8, 14, paste0("Lower than college (N=", 
                      format(edu.result.online$details$n.obs.1, big.mark = ","), 
                      ")"), pos = 4)
text(-1.8, 13, paste0("College or higher (N=", 
                      format(edu.result.online$details$n.obs.2, big.mark = ","), 
                      ")"), pos = 4)
text(-1.8, 12, "Difference b/w groups", pos = 4)
text(-1.95, 11, "Region", pos = 4, font = 2)
text(-1.8, 10, paste0("Luzon (N=", 
                      format(LV.result.online$details$n.obs.1, big.mark = ","), 
                      ")"), pos = 4)
text(-1.8, 9, paste0("Visayas (N=", 
                     format(LV.result.online$details$n.obs.2, big.mark = ","), 
                     ")"), pos = 4)
text(-1.8, 8, paste0("Mindanao (N=", 
                     format(VM.result.online$details$n.obs.2, big.mark = ","), 
                     ")"), pos = 4)
text(-1.8, 7, "Difference b/w L & V", pos = 4)
text(-1.8, 6, "Difference b/w L & M", pos = 4)
text(-1.8, 5, "Difference b/w V & M", pos = 4)
axis(1, seq(-0.5, 1, 0.25), lwd = 0.5, labels = NA)
mtext(sprintf("%4.2f", seq(-0.5, 1, 0.25)), 
      side = 1, line = 1, at = seq(-0.5, 1, 0.25))
mtext("Magnitude of SDB (or its difference)", side = 1, line = 2.5, at = 0.25)
dev.off()

#### estimated values ####
## gender
FtF.male.data <- subset(FtF.data, female == 0)
FtF.female.data <- subset(FtF.data, female == 1)
online.male.data <- subset(online.data, female == 0)
online.female.data <- subset(online.data, female == 1)

# DQ-based estimate (FtF, male)
DQ.est.FtF.male <- mean(FtF.male.data$DQ, na.rm = TRUE)
DQ.se.FtF.male <- sqrt(DQ.est.FtF.male * (1 - DQ.est.FtF.male) / 
                         sum(! is.na(FtF.male.data$DQ)))
round(DQ.est.FtF.male, 3)
round(DQ.est.FtF.male + qnorm(0.025) * DQ.se.FtF.male, 3)
round(DQ.est.FtF.male + qnorm(0.975) * DQ.se.FtF.male, 3)

# list-based estimate (FtF, male)
list.result.FtF.male <- lm_robust(list ~ treatment, data = FtF.male.data, 
                                  fixed_effects = ~ block)
list.est.FtF.male <- as.numeric(list.result.FtF.male$coefficients)
list.se.FtF.male <- as.numeric(list.result.FtF.male$std.error)
round(list.est.FtF.male, 3)
round(list.est.FtF.male + qnorm(0.025) * list.se.FtF.male, 3)
round(list.est.FtF.male + qnorm(0.975) * list.se.FtF.male, 3)

# SDB estimate (FtF, male)
SDB.est.FtF.male <- DQ.est.FtF.male - list.est.FtF.male
SDB.se.FtF.male <- sqrt(list.se.FtF.male ^ 2 + DQ.se.FtF.male ^ 2)
round(SDB.est.FtF.male, 3)
round(SDB.est.FtF.male + qnorm(0.025) * SDB.se.FtF.male, 3)
round(SDB.est.FtF.male + qnorm(0.975) * SDB.se.FtF.male, 3)
round((1 - pnorm(abs(SDB.est.FtF.male / SDB.se.FtF.male))) * 2, 3)

# DQ-based estimate (FtF, female)
DQ.est.FtF.female <- mean(FtF.female.data$DQ, na.rm = TRUE)
DQ.se.FtF.female <- sqrt(DQ.est.FtF.female * (1 - DQ.est.FtF.female) / 
                           sum(! is.na(FtF.female.data$DQ)))
round(DQ.est.FtF.female, 3)
round(DQ.est.FtF.female + qnorm(0.025) * DQ.se.FtF.female, 3)
round(DQ.est.FtF.female + qnorm(0.975) * DQ.se.FtF.female, 3)

# list-based estimate (FtF, female)
list.result.FtF.female <- lm_robust(list ~ treatment, data = FtF.female.data, 
                                    fixed_effects = ~ block)
list.est.FtF.female <- as.numeric(list.result.FtF.female$coefficients)
list.se.FtF.female <- as.numeric(list.result.FtF.female$std.error)
round(list.est.FtF.female, 3)
round(list.est.FtF.female + qnorm(0.025) * list.se.FtF.female, 3)
round(list.est.FtF.female + qnorm(0.975) * list.se.FtF.female, 3)

# SDB estimate (FtF, female)
SDB.est.FtF.female <- DQ.est.FtF.female - list.est.FtF.female
SDB.se.FtF.female <- sqrt(list.se.FtF.female ^ 2 + DQ.se.FtF.female ^ 2)
round(SDB.est.FtF.female, 3)
round(SDB.est.FtF.female + qnorm(0.025) * SDB.se.FtF.female, 3)
round(SDB.est.FtF.female + qnorm(0.975) * SDB.se.FtF.female, 3)
round((1 - pnorm(abs(SDB.est.FtF.female / SDB.se.FtF.female))) * 2, 3)

# DQ-based estimate (online, male)
DQ.est.online.male <- mean(online.male.data$DQ, na.rm = TRUE)
DQ.se.online.male <- sqrt(DQ.est.online.male * (1 - DQ.est.online.male) / 
                            sum(! is.na(online.male.data$DQ)))
round(DQ.est.online.male, 3)
round(DQ.est.online.male + qnorm(0.025) * DQ.se.online.male, 3)
round(DQ.est.online.male + qnorm(0.975) * DQ.se.online.male, 3)

# list-based estimate (online, male)
list.result.online.male <- lm_robust(list ~ treatment, data = online.male.data, 
                                     fixed_effects = ~ block)
list.est.online.male <- as.numeric(list.result.online.male$coefficients)
list.se.online.male <- as.numeric(list.result.online.male$std.error)
round(list.est.online.male, 3)
round(list.est.online.male + qnorm(0.025) * list.se.online.male, 3)
round(list.est.online.male + qnorm(0.975) * list.se.online.male, 3)

# SDB estimate (online, male)
SDB.est.online.male <- DQ.est.online.male - list.est.online.male
SDB.se.online.male <- sqrt(list.se.online.male ^ 2 + DQ.se.online.male ^ 2)
round(SDB.est.online.male, 3)
round(SDB.est.online.male + qnorm(0.025) * SDB.se.online.male, 3)
round(SDB.est.online.male + qnorm(0.975) * SDB.se.online.male, 3)
round((1 - pnorm(abs(SDB.est.online.male / SDB.se.online.male))) * 2, 3)

# DQ-based estimate (online, female)
DQ.est.online.female <- mean(online.female.data$DQ, na.rm = TRUE)
DQ.se.online.female <- sqrt(DQ.est.online.female * (1 - DQ.est.online.female) / 
                              sum(! is.na(online.female.data$DQ)))
round(DQ.est.online.female, 3)
round(DQ.est.online.female + qnorm(0.025) * DQ.se.online.female, 3)
round(DQ.est.online.female + qnorm(0.975) * DQ.se.online.female, 3)

# list-based estimate (online, female)
list.result.online.female <- lm_robust(list ~ treatment, data = online.female.data, 
                                       fixed_effects = ~ block)
list.est.online.female <- as.numeric(list.result.online.female$coefficients)
list.se.online.female <- as.numeric(list.result.online.female$std.error)
round(list.est.online.female, 3)
round(list.est.online.female + qnorm(0.025) * list.se.online.female, 3)
round(list.est.online.female + qnorm(0.975) * list.se.online.female, 3)

# SDB estimate (online, female)
SDB.est.online.female <- DQ.est.online.female - list.est.online.female
SDB.se.online.female <- sqrt(list.se.online.female ^ 2 + DQ.se.online.female ^ 2)
round(SDB.est.online.female, 3)
round(SDB.est.online.female + qnorm(0.025) * SDB.se.online.female, 3)
round(SDB.est.online.female + qnorm(0.975) * SDB.se.online.female, 3)
round((1 - pnorm(abs(SDB.est.online.female / SDB.se.online.female))) * 2, 3)

## education
FtF.lowedu.data <- subset(FtF.data, college == 0)
FtF.highedu.data <- subset(FtF.data, college == 1)
online.lowedu.data <- subset(online.data, college == 0)
online.highedu.data <- subset(online.data, college == 1)

# DQ-based estimate (FtF, low education)
DQ.est.FtF.lowedu <- mean(FtF.lowedu.data$DQ, na.rm = TRUE)
DQ.se.FtF.lowedu <- sqrt(DQ.est.FtF.lowedu * (1 - DQ.est.FtF.lowedu) / 
                           sum(! is.na(FtF.lowedu.data$DQ)))
round(DQ.est.FtF.lowedu, 3)
round(DQ.est.FtF.lowedu + qnorm(0.025) * DQ.se.FtF.lowedu, 3)
round(DQ.est.FtF.lowedu + qnorm(0.975) * DQ.se.FtF.lowedu, 3)

# list-based estimate (FtF, low education)
list.result.FtF.lowedu <- lm_robust(list ~ treatment, data = FtF.lowedu.data, 
                                    fixed_effects = ~ block)
list.est.FtF.lowedu <- as.numeric(list.result.FtF.lowedu$coefficients)
list.se.FtF.lowedu <- as.numeric(list.result.FtF.lowedu$std.error)
round(list.est.FtF.lowedu, 3)
round(list.est.FtF.lowedu + qnorm(0.025) * list.se.FtF.lowedu, 3)
round(list.est.FtF.lowedu + qnorm(0.975) * list.se.FtF.lowedu, 3)

# SDB estimate (FtF, low education)
SDB.est.FtF.lowedu <- DQ.est.FtF.lowedu - list.est.FtF.lowedu
SDB.se.FtF.lowedu <- sqrt(list.se.FtF.lowedu ^ 2 + DQ.se.FtF.lowedu ^ 2)
round(SDB.est.FtF.lowedu, 3)
round(SDB.est.FtF.lowedu + qnorm(0.025) * SDB.se.FtF.lowedu, 3)
round(SDB.est.FtF.lowedu + qnorm(0.975) * SDB.se.FtF.lowedu, 3)
round((1 - pnorm(abs(SDB.est.FtF.lowedu / SDB.se.FtF.lowedu))) * 2, 3)

# DQ-based estimate (FtF, high education)
DQ.est.FtF.highedu <- mean(FtF.highedu.data$DQ, na.rm = TRUE)
DQ.se.FtF.highedu <- sqrt(DQ.est.FtF.highedu * (1 - DQ.est.FtF.highedu) / 
                            sum(! is.na(FtF.highedu.data$DQ)))
round(DQ.est.FtF.highedu, 3)
round(DQ.est.FtF.highedu + qnorm(0.025) * DQ.se.FtF.highedu, 3)
round(DQ.est.FtF.highedu + qnorm(0.975) * DQ.se.FtF.highedu, 3)

# list-based estimate (FtF, high education)
list.result.FtF.highedu <- lm_robust(list ~ treatment, data = FtF.highedu.data, 
                                     fixed_effects = ~ block)
list.est.FtF.highedu <- as.numeric(list.result.FtF.highedu$coefficients)
list.se.FtF.highedu <- as.numeric(list.result.FtF.highedu$std.error)
round(list.est.FtF.highedu, 3)
round(list.est.FtF.highedu + qnorm(0.025) * list.se.FtF.highedu, 3)
round(list.est.FtF.highedu + qnorm(0.975) * list.se.FtF.highedu, 3)

# SDB estimate (FtF, high education)
SDB.est.FtF.highedu <- DQ.est.FtF.highedu - list.est.FtF.highedu
SDB.se.FtF.highedu <- sqrt(list.se.FtF.highedu ^ 2 + DQ.se.FtF.highedu ^ 2)
round(SDB.est.FtF.highedu, 3)
round(SDB.est.FtF.highedu + qnorm(0.025) * SDB.se.FtF.highedu, 3)
round(SDB.est.FtF.highedu + qnorm(0.975) * SDB.se.FtF.highedu, 3)
round((1 - pnorm(abs(SDB.est.FtF.highedu / SDB.se.FtF.highedu))) * 2, 3)

# DQ-based estimate (online, low education)
DQ.est.online.lowedu <- mean(online.lowedu.data$DQ, na.rm = TRUE)
DQ.se.online.lowedu <- sqrt(DQ.est.online.lowedu * (1 - DQ.est.online.lowedu) / 
                              sum(! is.na(online.lowedu.data$DQ)))
round(DQ.est.online.lowedu, 3)
round(DQ.est.online.lowedu + qnorm(0.025) * DQ.se.online.lowedu, 3)
round(DQ.est.online.lowedu + qnorm(0.975) * DQ.se.online.lowedu, 3)

# list-based estimate (online, low education)
list.result.online.lowedu <- lm_robust(list ~ treatment, data = online.lowedu.data, 
                                       fixed_effects = ~ block)
list.est.online.lowedu <- as.numeric(list.result.online.lowedu$coefficients)
list.se.online.lowedu <- as.numeric(list.result.online.lowedu$std.error)
round(list.est.online.lowedu, 3)
round(list.est.online.lowedu + qnorm(0.025) * list.se.online.lowedu, 3)
round(list.est.online.lowedu + qnorm(0.975) * list.se.online.lowedu, 3)

# SDB estimate (online, low education)
SDB.est.online.lowedu <- DQ.est.online.lowedu - list.est.online.lowedu
SDB.se.online.lowedu <- sqrt(list.se.online.lowedu ^ 2 + DQ.se.online.lowedu ^ 2)
round(SDB.est.online.lowedu, 3)
round(SDB.est.online.lowedu + qnorm(0.025) * SDB.se.online.lowedu, 3)
round(SDB.est.online.lowedu + qnorm(0.975) * SDB.se.online.lowedu, 3)
round((1 - pnorm(abs(SDB.est.online.lowedu / SDB.se.online.lowedu))) * 2, 3)

# DQ-based estimate (online, high education)
DQ.est.online.highedu <- mean(online.highedu.data$DQ, na.rm = TRUE)
DQ.se.online.highedu <- sqrt(DQ.est.online.highedu * (1 - DQ.est.online.highedu) / 
                               sum(! is.na(online.highedu.data$DQ)))
round(DQ.est.online.highedu, 3)
round(DQ.est.online.highedu + qnorm(0.025) * DQ.se.online.highedu, 3)
round(DQ.est.online.highedu + qnorm(0.975) * DQ.se.online.highedu, 3)

# list-based estimate (online, high education)
list.result.online.highedu <- lm_robust(list ~ treatment, data = online.highedu.data, 
                                        fixed_effects = ~ block)
list.est.online.highedu <- as.numeric(list.result.online.highedu$coefficients)
list.se.online.highedu <- as.numeric(list.result.online.highedu$std.error)
round(list.est.online.highedu, 3)
round(list.est.online.highedu + qnorm(0.025) * list.se.online.highedu, 3)
round(list.est.online.highedu + qnorm(0.975) * list.se.online.highedu, 3)

# SDB estimate (online, high education)
SDB.est.online.highedu <- DQ.est.online.highedu - list.est.online.highedu
SDB.se.online.highedu <- sqrt(list.se.online.highedu ^ 2 + DQ.se.online.highedu ^ 2)
round(SDB.est.online.highedu, 3)
round(SDB.est.online.highedu + qnorm(0.025) * SDB.se.online.highedu, 3)
round(SDB.est.online.highedu + qnorm(0.975) * SDB.se.online.highedu, 3)
round((1 - pnorm(abs(SDB.est.online.highedu / SDB.se.online.highedu))) * 2, 3)

## region
FtF.Luzon.data <- subset(FtF.data, Visayas == 0 & Mindanao == 0)
FtF.Visayas.data <- subset(FtF.data, Visayas == 1)
FtF.Mindanao.data <- subset(FtF.data, Mindanao == 1)
online.Luzon.data <- subset(online.data, Visayas == 0 & Mindanao == 0)
online.Visayas.data <- subset(online.data, Visayas == 1)
online.Mindanao.data <- subset(online.data, Mindanao == 1)

# DQ-based estimate (FtF, Luzon)
DQ.est.FtF.Luzon <- mean(FtF.Luzon.data$DQ, na.rm = TRUE)
DQ.se.FtF.Luzon <- sqrt(DQ.est.FtF.Luzon * (1 - DQ.est.FtF.Luzon) / sum(! is.na(FtF.Luzon.data$DQ)))
round(DQ.est.FtF.Luzon, 3)
round(DQ.est.FtF.Luzon + qnorm(0.025) * DQ.se.FtF.Luzon, 3)
round(DQ.est.FtF.Luzon + qnorm(0.975) * DQ.se.FtF.Luzon, 3)

# list-based estimate (FtF, Luzon)
list.result.FtF.Luzon <- lm_robust(list ~ treatment, data = FtF.Luzon.data,
                                   fixed_effects = ~ block)
list.est.FtF.Luzon <- as.numeric(list.result.FtF.Luzon$coefficients)
list.se.FtF.Luzon <- as.numeric(list.result.FtF.Luzon$std.error)
round(list.est.FtF.Luzon, 3)
round(list.est.FtF.Luzon + qnorm(0.025) * list.se.FtF.Luzon, 3)
round(list.est.FtF.Luzon + qnorm(0.975) * list.se.FtF.Luzon, 3)

# SDB estimate (FtF, Luzon)
SDB.est.FtF.Luzon <- DQ.est.FtF.Luzon - list.est.FtF.Luzon
SDB.se.FtF.Luzon <- sqrt(list.se.FtF.Luzon ^ 2 + DQ.se.FtF.Luzon ^ 2)
round(SDB.est.FtF.Luzon, 3)
round(SDB.est.FtF.Luzon + qnorm(0.025) * SDB.se.FtF.Luzon, 3)
round(SDB.est.FtF.Luzon + qnorm(0.975) * SDB.se.FtF.Luzon, 3)
round((1 - pnorm(abs(SDB.est.FtF.Luzon / SDB.se.FtF.Luzon))) * 2, 3)

# DQ-based estimate (FtF, Visayas)
DQ.est.FtF.Visayas <- mean(FtF.Visayas.data$DQ, na.rm = TRUE)
DQ.se.FtF.Visayas <- sqrt(DQ.est.FtF.Visayas * (1 - DQ.est.FtF.Visayas) / sum(! is.na(FtF.Visayas.data$DQ)))
round(DQ.est.FtF.Visayas, 3)
round(DQ.est.FtF.Visayas + qnorm(0.025) * DQ.se.FtF.Visayas, 3)
round(DQ.est.FtF.Visayas + qnorm(0.975) * DQ.se.FtF.Visayas, 3)

# list-based estimate (FtF, Visayas)
list.result.FtF.Visayas <- lm_robust(list ~ treatment, data = FtF.Visayas.data,
                                     fixed_effects = ~ block)
list.est.FtF.Visayas <- as.numeric(list.result.FtF.Visayas$coefficients)
list.se.FtF.Visayas <- as.numeric(list.result.FtF.Visayas$std.error)
round(list.est.FtF.Visayas, 3)
round(list.est.FtF.Visayas + qnorm(0.025) * list.se.FtF.Visayas, 3)
round(list.est.FtF.Visayas + qnorm(0.975) * list.se.FtF.Visayas, 3)

# SDB estimate (FtF, Visayas)
SDB.est.FtF.Visayas <- DQ.est.FtF.Visayas - list.est.FtF.Visayas
SDB.se.FtF.Visayas <- sqrt(list.se.FtF.Visayas ^ 2 + DQ.se.FtF.Visayas ^ 2)
round(SDB.est.FtF.Visayas, 3)
round(SDB.est.FtF.Visayas + qnorm(0.025) * SDB.se.FtF.Visayas, 3)
round(SDB.est.FtF.Visayas + qnorm(0.975) * SDB.se.FtF.Visayas, 3)
round((1 - pnorm(abs(SDB.est.FtF.Visayas / SDB.se.FtF.Visayas))) * 2, 3)

# DQ-based estimate (FtF, Mindanao)
DQ.est.FtF.Mindanao <- mean(FtF.Mindanao.data$DQ, na.rm = TRUE)
DQ.se.FtF.Mindanao <- sqrt(DQ.est.FtF.Mindanao * (1 - DQ.est.FtF.Mindanao) / sum(! is.na(FtF.Mindanao.data$DQ)))
round(DQ.est.FtF.Mindanao, 3)
round(DQ.est.FtF.Mindanao + qnorm(0.025) * DQ.se.FtF.Mindanao, 3)
round(DQ.est.FtF.Mindanao + qnorm(0.975) * DQ.se.FtF.Mindanao, 3)

# list-based estimate (FtF, Mindanao)
list.result.FtF.Mindanao <- lm_robust(list ~ treatment, data = FtF.Mindanao.data,
                                      fixed_effects = ~ block)
list.est.FtF.Mindanao <- as.numeric(list.result.FtF.Mindanao$coefficients)
list.se.FtF.Mindanao <- as.numeric(list.result.FtF.Mindanao$std.error)
round(list.est.FtF.Mindanao, 3)
round(list.est.FtF.Mindanao + qnorm(0.025) * list.se.FtF.Mindanao, 3)
round(list.est.FtF.Mindanao + qnorm(0.975) * list.se.FtF.Mindanao, 3)

# SDB estimate (FtF, Mindanao)
SDB.est.FtF.Mindanao <- DQ.est.FtF.Mindanao - list.est.FtF.Mindanao
SDB.se.FtF.Mindanao <- sqrt(list.se.FtF.Mindanao ^ 2 + DQ.se.FtF.Mindanao ^ 2)
round(SDB.est.FtF.Mindanao, 3)
round(SDB.est.FtF.Mindanao + qnorm(0.025) * SDB.se.FtF.Mindanao, 3)
round(SDB.est.FtF.Mindanao + qnorm(0.975) * SDB.se.FtF.Mindanao, 3)
round((1 - pnorm(abs(SDB.est.FtF.Mindanao / SDB.se.FtF.Mindanao))) * 2, 3)

# DQ-based estimate (online, Luzon)
DQ.est.online.Luzon <- mean(online.Luzon.data$DQ, na.rm = TRUE)
DQ.se.online.Luzon <- sqrt(DQ.est.online.Luzon * (1 - DQ.est.online.Luzon) / sum(! is.na(online.Luzon.data$DQ)))
round(DQ.est.online.Luzon, 3)
round(DQ.est.online.Luzon + qnorm(0.025) * DQ.se.online.Luzon, 3)
round(DQ.est.online.Luzon + qnorm(0.975) * DQ.se.online.Luzon, 3)

# list-based estimate (online, Luzon)
list.result.online.Luzon <- lm_robust(list ~ treatment, data = online.Luzon.data,
                                      fixed_effects = ~ block)
list.est.online.Luzon <- as.numeric(list.result.online.Luzon$coefficients)
list.se.online.Luzon <- as.numeric(list.result.online.Luzon$std.error)
round(list.est.online.Luzon, 3)
round(list.est.online.Luzon + qnorm(0.025) * list.se.online.Luzon, 3)
round(list.est.online.Luzon + qnorm(0.975) * list.se.online.Luzon, 3)

# SDB estimate (online, Luzon)
SDB.est.online.Luzon <- DQ.est.online.Luzon - list.est.online.Luzon
SDB.se.online.Luzon <- sqrt(list.se.online.Luzon ^ 2 + DQ.se.online.Luzon ^ 2)
round(SDB.est.online.Luzon, 3)
round(SDB.est.online.Luzon + qnorm(0.025) * SDB.se.online.Luzon, 3)
round(SDB.est.online.Luzon + qnorm(0.975) * SDB.se.online.Luzon, 3)
round((1 - pnorm(abs(SDB.est.online.Luzon / SDB.se.online.Luzon))) * 2, 3)

# DQ-based estimate (online, Visayas)
DQ.est.online.Visayas <- mean(online.Visayas.data$DQ, na.rm = TRUE)
DQ.se.online.Visayas <- sqrt(DQ.est.online.Visayas * (1 - DQ.est.online.Visayas) / sum(! is.na(online.Visayas.data$DQ)))
round(DQ.est.online.Visayas, 3)
round(DQ.est.online.Visayas + qnorm(0.025) * DQ.se.online.Visayas, 3)
round(DQ.est.online.Visayas + qnorm(0.975) * DQ.se.online.Visayas, 3)

# list-based estimate (online, Visayas)
list.result.online.Visayas <- lm_robust(list ~ treatment, data = online.Visayas.data,
                                        fixed_effects = ~ block)
list.est.online.Visayas <- as.numeric(list.result.online.Visayas$coefficients)
list.se.online.Visayas <- as.numeric(list.result.online.Visayas$std.error)
round(list.est.online.Visayas, 3)
round(list.est.online.Visayas + qnorm(0.025) * list.se.online.Visayas, 3)
round(list.est.online.Visayas + qnorm(0.975) * list.se.online.Visayas, 3)

# SDB estimate (online, Visayas)
SDB.est.online.Visayas <- DQ.est.online.Visayas - list.est.online.Visayas
SDB.se.online.Visayas <- sqrt(list.se.online.Visayas ^ 2 + DQ.se.online.Visayas ^ 2)
round(SDB.est.online.Visayas, 3)
round(SDB.est.online.Visayas + qnorm(0.025) * SDB.se.online.Visayas, 3)
round(SDB.est.online.Visayas + qnorm(0.975) * SDB.se.online.Visayas, 3)
round((1 - pnorm(abs(SDB.est.online.Visayas / SDB.se.online.Visayas))) * 2, 3)

# DQ-based estimate (online, Mindanao)
DQ.est.online.Mindanao <- mean(online.Mindanao.data$DQ, na.rm = TRUE)
DQ.se.online.Mindanao <- sqrt(DQ.est.online.Mindanao * (1 - DQ.est.online.Mindanao) / sum(! is.na(online.Mindanao.data$DQ)))
round(DQ.est.online.Mindanao, 3)
round(DQ.est.online.Mindanao + qnorm(0.025) * DQ.se.online.Mindanao, 3)
round(DQ.est.online.Mindanao + qnorm(0.975) * DQ.se.online.Mindanao, 3)

# list-based estimate (online, Mindanao)
list.result.online.Mindanao <- lm_robust(list ~ treatment, data = online.Mindanao.data,
                                         fixed_effects = ~ block)
list.est.online.Mindanao <- as.numeric(list.result.online.Mindanao$coefficients)
list.se.online.Mindanao <- as.numeric(list.result.online.Mindanao$std.error)
round(list.est.online.Mindanao, 3)
round(list.est.online.Mindanao + qnorm(0.025) * list.se.online.Mindanao, 3)
round(list.est.online.Mindanao + qnorm(0.975) * list.se.online.Mindanao, 3)

# SDB estimate (online, Mindanao)
SDB.est.online.Mindanao <- DQ.est.online.Mindanao - list.est.online.Mindanao
SDB.se.online.Mindanao <- sqrt(list.se.online.Mindanao ^ 2 + DQ.se.online.Mindanao ^ 2)
round(SDB.est.online.Mindanao, 3)
round(SDB.est.online.Mindanao + qnorm(0.025) * SDB.se.online.Mindanao, 3)
round(SDB.est.online.Mindanao + qnorm(0.975) * SDB.se.online.Mindanao, 3)
round((1 - pnorm(abs(SDB.est.online.Mindanao / SDB.se.online.Mindanao))) * 2, 3)

## socio-economic class
FtF.upper.data <- subset(FtF.data, lower == 0)
FtF.lower.data <- subset(FtF.data, lower == 1)

# DQ-based estimate (upper class)
DQ.est.FtF.upper <- mean(FtF.upper.data$DQ, na.rm = TRUE)
DQ.se.FtF.upper <- sqrt(DQ.est.FtF.upper * (1 - DQ.est.FtF.upper) / 
                          sum(! is.na(FtF.upper.data$DQ)))
round(DQ.est.FtF.upper, 3)
round(DQ.est.FtF.upper + qnorm(0.025) * DQ.se.FtF.upper, 3)
round(DQ.est.FtF.upper + qnorm(0.975) * DQ.se.FtF.upper, 3)

# list-based estimate (upper class)
list.result.FtF.upper <- lm_robust(list ~ treatment, data = FtF.upper.data, 
                                   fixed_effects = ~ block)
list.est.FtF.upper <- as.numeric(list.result.FtF.upper$coefficients)
list.se.FtF.upper <- as.numeric(list.result.FtF.upper$std.error)
round(list.est.FtF.upper, 3)
round(list.est.FtF.upper + qnorm(0.025) * list.se.FtF.upper, 3)
round(list.est.FtF.upper + qnorm(0.975) * list.se.FtF.upper, 3)

# SDB estimate (upper class)
SDB.est.FtF.upper <- DQ.est.FtF.upper - list.est.FtF.upper
SDB.se.FtF.upper <- sqrt(list.se.FtF.upper ^ 2 + DQ.se.FtF.upper ^ 2)
round(SDB.est.FtF.upper, 3)
round(SDB.est.FtF.upper + qnorm(0.025) * SDB.se.FtF.upper, 3)
round(SDB.est.FtF.upper + qnorm(0.975) * SDB.se.FtF.upper, 3)
round((1 - pnorm(abs(SDB.est.FtF.upper / SDB.se.FtF.upper))) * 2, 3)

# DQ-based estimate (lower class)
DQ.est.FtF.lower <- mean(FtF.lower.data$DQ, na.rm = TRUE)
DQ.se.FtF.lower <- sqrt(DQ.est.FtF.lower * (1 - DQ.est.FtF.lower) / 
                          sum(! is.na(FtF.lower.data$DQ)))
round(DQ.est.FtF.lower, 3)
round(DQ.est.FtF.lower + qnorm(0.025) * DQ.se.FtF.lower, 3)
round(DQ.est.FtF.lower + qnorm(0.975) * DQ.se.FtF.lower, 3)

# list-based estimate (lower class)
list.result.FtF.lower <- lm_robust(list ~ treatment, data = FtF.lower.data, 
                                   fixed_effects = ~ block)
list.est.FtF.lower <- as.numeric(list.result.FtF.lower$coefficients)
list.se.FtF.lower <- as.numeric(list.result.FtF.lower$std.error)
round(list.est.FtF.lower, 3)
round(list.est.FtF.lower + qnorm(0.025) * list.se.FtF.lower, 3)
round(list.est.FtF.lower + qnorm(0.975) * list.se.FtF.lower, 3)

# SDB estimate (lower class)
SDB.est.FtF.lower <- DQ.est.FtF.lower - list.est.FtF.lower
SDB.se.FtF.lower <- sqrt(list.se.FtF.lower ^ 2 + DQ.se.FtF.lower ^ 2)
round(SDB.est.FtF.lower, 3)
round(SDB.est.FtF.lower + qnorm(0.025) * SDB.se.FtF.lower, 3)
round(SDB.est.FtF.lower + qnorm(0.975) * SDB.se.FtF.lower, 3)
round((1 - pnorm(abs(SDB.est.FtF.lower / SDB.se.FtF.lower))) * 2, 3)

#### Figures A.9 and A.10 ####
cairo_pdf("Figure_A9.pdf", width = 6, height = 5.25, pointsize = 7)
layout(matrix(1:3, 3, 1))
par(mar = c(4, 0, 5, 1), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 6.4), ylim = c(0.5, 3.5), 
     main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = c(seq(0, 1, 0.25), seq(1.4, 2.9, 0.25), 
             seq(3.5, 4.5, 0.25), seq(4.9, 6.4, 0.25)), col = "gray")
segments(DQ.est.FtF.male + qnorm(0.025) * DQ.se.FtF.male, 3.1, 
         DQ.est.FtF.male + qnorm(0.975) * DQ.se.FtF.male, 3.1)
points(DQ.est.FtF.male, 3.1, pch = 15)
text(DQ.est.FtF.male, 3.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.FtF.male + qnorm(0.025) * list.se.FtF.male, 2.9, 
         list.est.FtF.male + qnorm(0.975) * list.se.FtF.male, 2.9)
points(list.est.FtF.male, 2.9, pch = 19)
text(list.est.FtF.male, 2.9, "List", pos = 1, cex = 0.8)
segments(DQ.est.FtF.female + qnorm(0.025) * DQ.se.FtF.female, 2.1, 
         DQ.est.FtF.female + qnorm(0.975) * DQ.se.FtF.female, 2.1)
points(DQ.est.FtF.female, 2.1, pch = 15)
text(DQ.est.FtF.female, 2.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.FtF.female + qnorm(0.025) * list.se.FtF.female, 1.9, 
         list.est.FtF.female + qnorm(0.975) * list.se.FtF.female, 1.9)
points(list.est.FtF.female, 1.9, pch = 19)
text(list.est.FtF.female, 1.9, "List", pos = 1, cex = 0.8)
segments(SDB.est.FtF.male + qnorm(0.025) * SDB.se.FtF.male + 1.9, 3, 
         SDB.est.FtF.male + qnorm(0.975) * SDB.se.FtF.male + 1.9, 3)
points(SDB.est.FtF.male + 1.9, 3, pch = 4)
text(SDB.est.FtF.male + 1.9, 3, sprintf("%5.3f", SDB.est.FtF.male), 
     cex = 0.8, pos = 3)
segments(SDB.est.FtF.female + qnorm(0.025) * SDB.se.FtF.female + 1.9, 2, 
         SDB.est.FtF.female + qnorm(0.975) * SDB.se.FtF.female + 1.9, 2)
points(SDB.est.FtF.female + 1.9, 2, pch = 4)
text(SDB.est.FtF.female + 1.9, 2, sprintf("%5.3f", SDB.est.FtF.female), 
     cex = 0.8, pos = 3)
segments(gender.result.FtF$details$SDB.diff.est + qnorm(0.025) * 
           gender.result.FtF$details$SDB.diff.se + 1.9, 1, 
         gender.result.FtF$details$SDB.diff.est + qnorm(0.975) * 
           gender.result.FtF$details$SDB.diff.se + 1.9, 1)
points(gender.result.FtF$details$SDB.diff.est + 1.9, 1, pch = 4)
text(gender.result.FtF$details$SDB.diff.est + 1.9, 1, 
     sprintf("%5.3f", gender.result.FtF$details$SDB.diff.est), 
     cex = 0.8, pos = 3)
segments(DQ.est.online.male + qnorm(0.025) * DQ.se.online.male + 3.5, 3.1, 
         DQ.est.online.male + qnorm(0.975) * DQ.se.online.male + 3.5, 3.1)
points(DQ.est.online.male + 3.5, 3.1, pch = 15)
text(DQ.est.online.male + 3.5, 3.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.online.male + qnorm(0.025) * list.se.online.male + 3.5, 2.9, 
         list.est.online.male + qnorm(0.975) * list.se.online.male + 3.5, 2.9)
points(list.est.online.male + 3.5, 2.9, pch = 19)
text(list.est.online.male + 3.5, 2.9, "List", pos = 1, cex = 0.8)
segments(DQ.est.online.female + qnorm(0.025) * DQ.se.online.female + 3.5, 2.1, 
         DQ.est.online.female + qnorm(0.975) * DQ.se.online.female + 3.5, 2.1)
points(DQ.est.online.female + 3.5, 2.1, pch = 15)
text(DQ.est.online.female + 3.5, 2.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.online.female + qnorm(0.025) * list.se.online.female + 3.5, 1.9, 
         list.est.online.female + qnorm(0.975) * list.se.online.female + 3.5, 1.9)
points(list.est.online.female + 3.5, 1.9, pch = 19)
text(list.est.online.female + 3.5, 1.9, "List", pos = 1, cex = 0.8)
segments(SDB.est.online.male + qnorm(0.025) * SDB.se.online.male + 5.4, 3, 
         SDB.est.online.male + qnorm(0.975) * SDB.se.online.male + 5.4, 3)
points(SDB.est.online.male + 5.4, 3, pch = 4)
text(SDB.est.online.male + 5.4, 3, sprintf("%5.3f", SDB.est.online.male), 
     cex = 0.8, pos = 3)
segments(SDB.est.online.female + qnorm(0.025) * SDB.se.online.female + 5.4, 2, 
         SDB.est.online.female + qnorm(0.975) * SDB.se.online.female + 5.4, 2)
points(SDB.est.online.female + 5.4, 2, pch = 4)
text(SDB.est.online.female + 5.4, 2, sprintf("%5.3f", SDB.est.online.female), 
     cex = 0.8, pos = 3)
segments(gender.result.online$details$SDB.diff.est + qnorm(0.025) * 
           gender.result.online$details$SDB.diff.se + 5.4, 1, 
         gender.result.online$details$SDB.diff.est + qnorm(0.975) * 
           gender.result.online$details$SDB.diff.se + 5.4, 1)
points(gender.result.online$details$SDB.diff.est + 5.4, 1, pch = 4)
text(gender.result.online$details$SDB.diff.est + 5.4, 1, 
     sprintf("%5.3f", gender.result.online$details$SDB.diff.est), 
     cex = 0.8, pos = 3)
text(-1.2, 3, "Male", pos = 4, cex = 1.2)
text(-1.2, 2, "Female", pos = 4, cex = 1.2)
text(-1.2, 1, "Difference b/w groups", pos = 4, cex = 1.2)
axis(1, seq(0, 1, 0.25), labels = NA, lwd = 0.5)
mtext(sprintf("%4.2f", seq(0, 1, 0.25)), side = 1, line = 1, 
      at = seq(0, 1, 0.25), cex = 0.7)
axis(1, seq(1.4, 2.9, 0.25), labels = NA, lwd = 0.5)
mtext(sprintf("%4.2f", seq(-0.5, 1, 0.25)), side = 1, line = 1, 
      at = seq(1.4, 2.9, 0.25), cex = 0.7)
axis(1, seq(3.5, 4.5, 0.25), labels = NA, lwd = 0.5)
mtext(sprintf("%4.2f", seq(0, 1, 0.25)), side = 1, line = 1, 
      at = seq(3.5, 4.5, 0.25), cex = 0.7)
axis(1, seq(4.9, 6.4, 0.25), labels = NA, lwd = 0.5)
mtext(sprintf("%4.2f", seq(-0.5, 1, 0.25)), side = 1, line = 1, 
      at = seq(4.9, 6.4, 0.25), cex = 0.7)
mtext("Percent", side = 1, line = 2.5, at = 0.5, cex = 0.8)
mtext("Percentage points", side = 1, line = 2.5, at = 2.15, cex = 0.8)
mtext("Percent", side = 1, line = 2.5, at = 4, cex = 0.8)
mtext("Percentage points", side = 1, line = 2.5, at = 5.65, cex = 0.8)
mtext("Gender", line = 2.5, at = 3.2, font = 2, cex = 1.3)
mtext("Face-to-face", line = 1.5, at = 1.45, font = 2, cex = 1.1)
mtext("Online", line = 1.5, at = 4.95, font = 2, cex = 1.1)
mtext("Approval rate", line = 0.5, at = 0.5, font = 2, cex = 0.9)
mtext("SDB", line = 0.5, at = 2.15, font = 2, cex = 0.9)
mtext("Approval rate", line = 0.5, at = 4, font = 2, cex = 0.9)
mtext("SDB", line = 0.5, at = 5.65, font = 2, cex = 0.9)
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 6.4), ylim = c(0.5, 3.5), 
     main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = c(seq(0, 1, 0.25), seq(1.4, 2.9, 0.25), 
             seq(3.5, 4.5, 0.25), seq(4.9, 6.4, 0.25)), col = "gray")
segments(DQ.est.FtF.lowedu + qnorm(0.025) * DQ.se.FtF.lowedu, 3.1, 
         DQ.est.FtF.lowedu + qnorm(0.975) * DQ.se.FtF.lowedu, 3.1)
points(DQ.est.FtF.lowedu, 3.1, pch = 15)
text(DQ.est.FtF.lowedu, 3.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.FtF.lowedu + qnorm(0.025) * list.se.FtF.lowedu, 2.9, 
         list.est.FtF.lowedu + qnorm(0.975) * list.se.FtF.lowedu, 2.9)
points(list.est.FtF.lowedu, 2.9, pch = 19)
text(list.est.FtF.lowedu, 2.9, "List", pos = 1, cex = 0.8)
segments(DQ.est.FtF.highedu + qnorm(0.025) * DQ.se.FtF.highedu, 2.1, 
         DQ.est.FtF.highedu + qnorm(0.975) * DQ.se.FtF.highedu, 2.1)
points(DQ.est.FtF.highedu, 2.1, pch = 15)
text(DQ.est.FtF.highedu, 2.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.FtF.highedu + qnorm(0.025) * list.se.FtF.highedu, 1.9, 
         list.est.FtF.highedu + qnorm(0.975) * list.se.FtF.highedu, 1.9)
points(list.est.FtF.highedu, 1.9, pch = 19)
text(list.est.FtF.highedu, 1.9, "List", pos = 1, cex = 0.8)
segments(SDB.est.FtF.lowedu + qnorm(0.025) * SDB.se.FtF.lowedu + 1.9, 3, 
         SDB.est.FtF.lowedu + qnorm(0.975) * SDB.se.FtF.lowedu + 1.9, 3)
points(SDB.est.FtF.lowedu + 1.9, 3, pch = 4)
text(SDB.est.FtF.lowedu + 1.9, 3, sprintf("%5.3f", SDB.est.FtF.lowedu), 
     cex = 0.8, pos = 3)
segments(SDB.est.FtF.highedu + qnorm(0.025) * SDB.se.FtF.highedu + 1.9, 2, 
         SDB.est.FtF.highedu + qnorm(0.975) * SDB.se.FtF.highedu + 1.9, 2)
points(SDB.est.FtF.highedu + 1.9, 2, pch = 4)
text(SDB.est.FtF.highedu + 1.9, 2, sprintf("%5.3f", SDB.est.FtF.highedu), 
     cex = 0.8, pos = 3)
segments(edu.result.FtF$details$SDB.diff.est + qnorm(0.025) * 
           edu.result.FtF$details$SDB.diff.se + 1.9, 1, 
         edu.result.FtF$details$SDB.diff.est + qnorm(0.975) * 
           edu.result.FtF$details$SDB.diff.se + 1.9, 1)
points(edu.result.FtF$details$SDB.diff.est + 1.9, 1, pch = 4)
text(edu.result.FtF$details$SDB.diff.est + 1.9, 1, 
     sprintf("%5.3f", edu.result.FtF$details$SDB.diff.est), 
     cex = 0.8, pos = 3)
segments(DQ.est.online.lowedu + qnorm(0.025) * DQ.se.online.lowedu + 3.5, 3.1, 
         DQ.est.online.lowedu + qnorm(0.975) * DQ.se.online.lowedu + 3.5, 3.1)
points(DQ.est.online.lowedu + 3.5, 3.1, pch = 15)
text(DQ.est.online.lowedu + 3.5, 3.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.online.lowedu + qnorm(0.025) * list.se.online.lowedu + 3.5, 2.9, 
         list.est.online.lowedu + qnorm(0.975) * list.se.online.lowedu + 3.5, 2.9)
points(list.est.online.lowedu + 3.5, 2.9, pch = 19)
text(list.est.online.lowedu + 3.5, 2.9, "List", pos = 1, cex = 0.8)
segments(DQ.est.online.highedu + qnorm(0.025) * DQ.se.online.highedu + 3.5, 2.1, 
         DQ.est.online.highedu + qnorm(0.975) * DQ.se.online.highedu + 3.5, 2.1)
points(DQ.est.online.highedu + 3.5, 2.1, pch = 15)
text(DQ.est.online.highedu + 3.5, 2.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.online.highedu + qnorm(0.025) * list.se.online.highedu + 3.5, 1.9, 
         list.est.online.highedu + qnorm(0.975) * list.se.online.highedu + 3.5, 1.9)
points(list.est.online.highedu + 3.5, 1.9, pch = 19)
text(list.est.online.highedu + 3.5, 1.9, "List", pos = 1, cex = 0.8)
segments(SDB.est.online.lowedu + qnorm(0.025) * SDB.se.online.lowedu + 5.4, 3, 
         SDB.est.online.lowedu + qnorm(0.975) * SDB.se.online.lowedu + 5.4, 3)
points(SDB.est.online.lowedu + 5.4, 3, pch = 4)
text(SDB.est.online.lowedu + 5.4, 3, sprintf("%5.3f", SDB.est.online.lowedu), 
     cex = 0.8, pos = 3)
segments(SDB.est.online.highedu + qnorm(0.025) * SDB.se.online.highedu + 5.4, 2, 
         SDB.est.online.highedu + qnorm(0.975) * SDB.se.online.highedu + 5.4, 2)
points(SDB.est.online.highedu + 5.4, 2, pch = 4)
text(SDB.est.online.highedu + 5.4, 2, sprintf("%5.3f", SDB.est.online.highedu), 
     cex = 0.8, pos = 3)
segments(edu.result.online$details$SDB.diff.est + qnorm(0.025) * 
           edu.result.online$details$SDB.diff.se + 5.4, 1, 
         edu.result.online$details$SDB.diff.est + qnorm(0.975) * 
           edu.result.online$details$SDB.diff.se + 5.4, 1)
points(edu.result.online$details$SDB.diff.est + 5.4, 1, pch = 4)
text(edu.result.online$details$SDB.diff.est + 5.4, 1, 
     sprintf("%5.3f", edu.result.online$details$SDB.diff.est), 
     cex = 0.8, pos = 3)
text(-1.2, 3, "Lower than college", pos = 4, cex = 1.2)
text(-1.2, 2, "College or higher", pos = 4, cex = 1.2)
text(-1.2, 1, "Difference b/w groups", pos = 4, cex = 1.2)
axis(1, seq(0, 1, 0.25), labels = NA, lwd = 0.5)
mtext(sprintf("%4.2f", seq(0, 1, 0.25)), side = 1, line = 1, 
      at = seq(0, 1, 0.25), cex = 0.7)
axis(1, seq(1.4, 2.9, 0.25), labels = NA, lwd = 0.5)
mtext(sprintf("%4.2f", seq(-0.5, 1, 0.25)), side = 1, line = 1, 
      at = seq(1.4, 2.9, 0.25), cex = 0.7)
axis(1, seq(3.5, 4.5, 0.25), labels = NA, lwd = 0.5)
mtext(sprintf("%4.2f", seq(0, 1, 0.25)), side = 1, line = 1, 
      at = seq(3.5, 4.5, 0.25), cex = 0.7)
axis(1, seq(4.9, 6.4, 0.25), labels = NA, lwd = 0.5)
mtext(sprintf("%4.2f", seq(-0.5, 1, 0.25)), side = 1, line = 1, 
      at = seq(4.9, 6.4, 0.25), cex = 0.7)
mtext("Percent", side = 1, line = 2.5, at = 0.5, cex = 0.8)
mtext("Percentage points", side = 1, line = 2.5, at = 2.15, cex = 0.8)
mtext("Percent", side = 1, line = 2.5, at = 4, cex = 0.8)
mtext("Percentage points", side = 1, line = 2.5, at = 5.65, cex = 0.8)
mtext("Education", line = 2.5, at = 3.2, font = 2, cex = 1.3)
mtext("Face-to-face", line = 1.5, at = 1.45, font = 2, cex = 1.1)
mtext("Online", line = 1.5, at = 4.95, font = 2, cex = 1.1)
mtext("Approval rate", line = 0.5, at = 0.5, font = 2, cex = 0.9)
mtext("SDB", line = 0.5, at = 2.15, font = 2, cex = 0.9)
mtext("Approval rate", line = 0.5, at = 4, font = 2, cex = 0.9)
mtext("SDB", line = 0.5, at = 5.65, font = 2, cex = 0.9)
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 6.4), ylim = c(0.5, 3.5), 
     main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = c(seq(0, 1, 0.25), seq(1.4, 2.9, 0.25)), col = "gray")
segments(DQ.est.FtF.upper + qnorm(0.025) * DQ.se.FtF.upper, 3.1, 
         DQ.est.FtF.upper + qnorm(0.975) * DQ.se.FtF.upper, 3.1)
points(DQ.est.FtF.upper, 3.1, pch = 15)
text(DQ.est.FtF.upper, 3.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.FtF.upper + qnorm(0.025) * list.se.FtF.upper, 2.9, 
         list.est.FtF.upper + qnorm(0.975) * list.se.FtF.upper, 2.9)
points(list.est.FtF.upper, 2.9, pch = 19)
text(list.est.FtF.upper, 2.9, "List", pos = 1, cex = 0.8)
segments(DQ.est.FtF.lower + qnorm(0.025) * DQ.se.FtF.lower, 2.1, 
         DQ.est.FtF.lower + qnorm(0.975) * DQ.se.FtF.lower, 2.1)
points(DQ.est.FtF.lower, 2.1, pch = 15)
text(DQ.est.FtF.lower, 2.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.FtF.lower + qnorm(0.025) * list.se.FtF.lower, 1.9, 
         list.est.FtF.lower + qnorm(0.975) * list.se.FtF.lower, 1.9)
points(list.est.FtF.lower, 1.9, pch = 19)
text(list.est.FtF.lower, 1.9, "List", pos = 1, cex = 0.8)
segments(SDB.est.FtF.upper + qnorm(0.025) * SDB.se.FtF.upper + 1.9, 3, 
         SDB.est.FtF.upper + qnorm(0.975) * SDB.se.FtF.upper + 1.9, 3)
points(SDB.est.FtF.upper + 1.9, 3, pch = 4)
text(SDB.est.FtF.upper + 1.9, 3, sprintf("%5.3f", SDB.est.FtF.upper), 
     cex = 0.8, pos = 3)
segments(SDB.est.FtF.lower + qnorm(0.025) * SDB.se.FtF.lower + 1.9, 2, 
         SDB.est.FtF.lower + qnorm(0.975) * SDB.se.FtF.lower + 1.9, 2)
points(SDB.est.FtF.lower + 1.9, 2, pch = 4)
text(SDB.est.FtF.lower + 1.9, 2, sprintf("%5.3f", SDB.est.FtF.lower), 
     cex = 0.8, pos = 3)
segments(class.result.FtF$details$SDB.diff.est + qnorm(0.025) * 
           class.result.FtF$details$SDB.diff.se + 1.9, 1, 
         class.result.FtF$details$SDB.diff.est + qnorm(0.975) * 
           class.result.FtF$details$SDB.diff.se + 1.9, 1)
points(class.result.FtF$details$SDB.diff.est + 1.9, 1, pch = 4)
text(class.result.FtF$details$SDB.diff.est + 1.9, 1, 
     sprintf("%5.3f", class.result.FtF$details$SDB.diff.est), 
     cex = 0.8, pos = 3)
text(-1.2, 3, "Upper class", pos = 4, cex = 1.2)
text(-1.2, 2, "Lower class", pos = 4, cex = 1.2)
text(-1.2, 1, "Difference b/w groups", pos = 4, cex = 1.2)
axis(1, seq(0, 1, 0.25), labels = NA, lwd = 0.5)
mtext(sprintf("%4.2f", seq(0, 1, 0.25)), side = 1, line = 1, 
      at = seq(0, 1, 0.25), cex = 0.7)
axis(1, seq(1.4, 2.9, 0.25), labels = NA, lwd = 0.5)
mtext(sprintf("%4.2f", seq(-0.5, 1, 0.25)), side = 1, line = 1, 
      at = seq(1.4, 2.9, 0.25), cex = 0.7)
mtext("Percent", side = 1, line = 2.5, at = 0.5, cex = 0.8)
mtext("Percentage points", side = 1, line = 2.5, at = 2.15, cex = 0.8)
mtext("Socio-economic class", line = 2.5, at = 3.2, font = 2, cex = 1.3)
mtext("Face-to-face", line = 1.5, at = 1.45, font = 2, cex = 1.1)
mtext("Approval rate", line = 0.5, at = 0.5, font = 2, cex = 0.9)
mtext("SDB", line = 0.5, at = 2.15, font = 2, cex = 0.9)
dev.off()

cairo_pdf("Figure_A10.pdf", width = 6, height = 3.2, pointsize = 6)
par(mar = c(4, 0, 5, 1), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-1.6, 6.4), ylim = c(0.5, 6.5), 
     main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = c(seq(0, 1, 0.25), seq(1.4, 2.9, 0.25), 
             seq(3.5, 4.5, 0.25), seq(4.9, 6.4, 0.25)), col = "gray")
segments(DQ.est.FtF.Luzon + qnorm(0.025) * DQ.se.FtF.Luzon, 6.1, 
         DQ.est.FtF.Luzon + qnorm(0.975) * DQ.se.FtF.Luzon, 6.1)
points(DQ.est.FtF.Luzon, 6.1, pch = 15)
text(DQ.est.FtF.Luzon, 6.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.FtF.Luzon + qnorm(0.025) * list.se.FtF.Luzon, 5.9, 
         list.est.FtF.Luzon + qnorm(0.975) * list.se.FtF.Luzon, 5.9)
points(list.est.FtF.Luzon, 5.9, pch = 19)
text(list.est.FtF.Luzon, 5.9, "List", pos = 1, cex = 0.8)
segments(DQ.est.FtF.Visayas + qnorm(0.025) * DQ.se.FtF.Visayas, 5.1, 
         DQ.est.FtF.Visayas + qnorm(0.975) * DQ.se.FtF.Visayas, 5.1)
points(DQ.est.FtF.Visayas, 5.1, pch = 15)
text(DQ.est.FtF.Visayas, 5.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.FtF.Visayas + qnorm(0.025) * list.se.FtF.Visayas, 4.9, 
         list.est.FtF.Visayas + qnorm(0.975) * list.se.FtF.Visayas, 4.9)
points(list.est.FtF.Visayas, 4.9, pch = 19)
text(list.est.FtF.Visayas, 4.9, "List", pos = 1, cex = 0.8)
segments(DQ.est.FtF.Mindanao + qnorm(0.025) * DQ.se.FtF.Mindanao, 4.1, 
         DQ.est.FtF.Mindanao + qnorm(0.975) * DQ.se.FtF.Mindanao, 4.1)
points(DQ.est.FtF.Mindanao, 4.1, pch = 15)
text(DQ.est.FtF.Mindanao, 4.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.FtF.Mindanao + qnorm(0.025) * list.se.FtF.Mindanao, 3.9, 
         list.est.FtF.Mindanao + qnorm(0.975) * list.se.FtF.Mindanao, 3.9)
points(list.est.FtF.Mindanao, 3.9, pch = 19)
text(list.est.FtF.Mindanao, 3.9, "List", pos = 1, cex = 0.8)
segments(SDB.est.FtF.Luzon + qnorm(0.025) * SDB.se.FtF.Luzon + 1.9, 6, 
         SDB.est.FtF.Luzon + qnorm(0.975) * SDB.se.FtF.Luzon + 1.9, 6)
points(SDB.est.FtF.Luzon + 1.9, 6, pch = 4)
text(SDB.est.FtF.Luzon + 1.9, 6, sprintf("%5.3f", SDB.est.FtF.Luzon), 
     cex = 0.8, pos = 3)
segments(SDB.est.FtF.Visayas + qnorm(0.025) * SDB.se.FtF.Visayas + 1.9, 5, 
         SDB.est.FtF.Visayas + qnorm(0.975) * SDB.se.FtF.Visayas + 1.9, 5)
points(SDB.est.FtF.Visayas + 1.9, 5, pch = 4)
text(SDB.est.FtF.Visayas + 1.9, 5, sprintf("%5.3f", SDB.est.FtF.Visayas), 
     cex = 0.8, pos = 3)
segments(SDB.est.FtF.Mindanao + qnorm(0.025) * SDB.se.FtF.Mindanao + 1.9, 4, 
         SDB.est.FtF.Mindanao + qnorm(0.975) * SDB.se.FtF.Mindanao + 1.9, 4)
points(SDB.est.FtF.Mindanao + 1.9, 4, pch = 4)
text(SDB.est.FtF.Mindanao + 1.9, 4, sprintf("%5.3f", SDB.est.FtF.Mindanao), 
     cex = 0.8, pos = 3)
segments(LV.result.FtF$details$SDB.diff.est + qnorm(0.025) * 
           LV.result.FtF$details$SDB.diff.se + 1.9, 3, 
         LV.result.FtF$details$SDB.diff.est + qnorm(0.975) * 
           LV.result.FtF$details$SDB.diff.se + 1.9, 3)
points(LV.result.FtF$details$SDB.diff.est + 1.9, 3, pch = 4)
text(LV.result.FtF$details$SDB.diff.est + 1.9, 3, 
     sprintf("%5.3f", LV.result.FtF$details$SDB.diff.est), 
     cex = 0.8, pos = 3)
segments(LM.result.FtF$details$SDB.diff.est + qnorm(0.025) * 
           LM.result.FtF$details$SDB.diff.se + 1.9, 2, 
         LM.result.FtF$details$SDB.diff.est + qnorm(0.975) * 
           LM.result.FtF$details$SDB.diff.se + 1.9, 2)
points(LM.result.FtF$details$SDB.diff.est + 1.9, 2, pch = 4)
text(LM.result.FtF$details$SDB.diff.est + 1.9, 2, 
     sprintf("%5.3f", LM.result.FtF$details$SDB.diff.est), 
     cex = 0.8, pos = 3)
segments(VM.result.FtF$details$SDB.diff.est + qnorm(0.025) * 
           VM.result.FtF$details$SDB.diff.se + 1.9, 1, 
         VM.result.FtF$details$SDB.diff.est + qnorm(0.975) * 
           VM.result.FtF$details$SDB.diff.se + 1.9, 1)
points(VM.result.FtF$details$SDB.diff.est + 1.9, 1, pch = 4)
text(VM.result.FtF$details$SDB.diff.est + 1.9, 1, 
     sprintf("%5.3f", VM.result.FtF$details$SDB.diff.est), 
     cex = 0.8, pos = 3)
segments(DQ.est.online.Luzon + qnorm(0.025) * DQ.se.online.Luzon + 3.5, 6.1, 
         DQ.est.online.Luzon + qnorm(0.975) * DQ.se.online.Luzon + 3.5, 6.1)
points(DQ.est.online.Luzon + 3.5, 6.1, pch = 15)
text(DQ.est.online.Luzon + 3.5, 6.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.online.Luzon + qnorm(0.025) * list.se.online.Luzon + 3.5, 5.9, 
         list.est.online.Luzon + qnorm(0.975) * list.se.online.Luzon + 3.5, 5.9)
points(list.est.online.Luzon + 3.5, 5.9, pch = 19)
text(list.est.online.Luzon + 3.5, 5.9, "List", pos = 1, cex = 0.8)
segments(DQ.est.online.Visayas + qnorm(0.025) * DQ.se.online.Visayas + 3.5, 5.1, 
         DQ.est.online.Visayas + qnorm(0.975) * DQ.se.online.Visayas + 3.5, 5.1)
points(DQ.est.online.Visayas + 3.5, 5.1, pch = 15)
text(DQ.est.online.Visayas + 3.5, 5.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.online.Visayas + qnorm(0.025) * list.se.online.Visayas + 3.5, 4.9, 
         list.est.online.Visayas + qnorm(0.975) * list.se.online.Visayas + 3.5, 4.9)
points(list.est.online.Visayas + 3.5, 4.9, pch = 19)
text(list.est.online.Visayas + 3.5, 4.9, "List", pos = 1, cex = 0.8)
segments(DQ.est.online.Mindanao + qnorm(0.025) * DQ.se.online.Mindanao + 3.5, 4.1, 
         DQ.est.online.Mindanao + qnorm(0.975) * DQ.se.online.Mindanao + 3.5, 4.1)
points(DQ.est.online.Mindanao + 3.5, 4.1, pch = 15)
text(DQ.est.online.Mindanao + 3.5, 4.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.online.Mindanao + qnorm(0.025) * list.se.online.Mindanao + 3.5, 3.9, 
         list.est.online.Mindanao + qnorm(0.975) * list.se.online.Mindanao + 3.5, 3.9)
points(list.est.online.Mindanao + 3.5, 3.9, pch = 19)
text(list.est.online.Mindanao + 3.5, 3.9, "List", pos = 1, cex = 0.8)
segments(SDB.est.online.Luzon + qnorm(0.025) * SDB.se.online.Luzon + 5.4, 6, 
         SDB.est.online.Luzon + qnorm(0.975) * SDB.se.online.Luzon + 5.4, 6)
points(SDB.est.online.Luzon + 5.4, 6, pch = 4)
text(SDB.est.online.Luzon + 5.4, 6, sprintf("%5.3f", SDB.est.online.Luzon), 
     cex = 0.8, pos = 3)
segments(SDB.est.online.Visayas + qnorm(0.025) * SDB.se.online.Visayas + 5.4, 5, 
         SDB.est.online.Visayas + qnorm(0.975) * SDB.se.online.Visayas + 5.4, 5)
points(SDB.est.online.Visayas + 5.4, 5, pch = 4)
text(SDB.est.online.Visayas + 5.4, 5, sprintf("%5.3f", SDB.est.online.Visayas), 
     cex = 0.8, pos = 3)
segments(SDB.est.online.Mindanao + qnorm(0.025) * SDB.se.online.Mindanao + 5.4, 4, 
         SDB.est.online.Mindanao + qnorm(0.975) * SDB.se.online.Mindanao + 5.4, 4)
points(SDB.est.online.Mindanao + 5.4, 4, pch = 4)
text(SDB.est.online.Mindanao + 5.4, 4, sprintf("%5.3f", SDB.est.online.Mindanao), 
     cex = 0.8, pos = 3)
segments(LV.result.online$details$SDB.diff.est + qnorm(0.025) * 
           LV.result.online$details$SDB.diff.se + 5.4, 3, 
         LV.result.online$details$SDB.diff.est + qnorm(0.975) * 
           LV.result.online$details$SDB.diff.se + 5.4, 3)
points(LV.result.online$details$SDB.diff.est + 5.4, 3, pch = 4)
text(LV.result.online$details$SDB.diff.est + 5.4, 3, 
     sprintf("%5.3f", LV.result.online$details$SDB.diff.est), 
     cex = 0.8, pos = 3)
segments(LM.result.online$details$SDB.diff.est + qnorm(0.025) * 
           LM.result.online$details$SDB.diff.se + 5.4, 2, 
         LM.result.online$details$SDB.diff.est + qnorm(0.975) * 
           LM.result.online$details$SDB.diff.se + 5.4, 2)
points(LM.result.online$details$SDB.diff.est + 5.4, 2, pch = 4)
text(LM.result.online$details$SDB.diff.est + 5.4, 2, 
     sprintf("%5.3f", LM.result.online$details$SDB.diff.est), 
     cex = 0.8, pos = 3)
segments(VM.result.online$details$SDB.diff.est + qnorm(0.025) * 
           VM.result.online$details$SDB.diff.se + 5.4, 1, 
         VM.result.online$details$SDB.diff.est + qnorm(0.975) * 
           VM.result.online$details$SDB.diff.se + 5.4, 1)
points(VM.result.online$details$SDB.diff.est + 5.4, 1, pch = 4)
text(VM.result.online$details$SDB.diff.est + 5.4, 1, 
     sprintf("%5.3f", VM.result.online$details$SDB.diff.est), 
     cex = 0.8, pos = 3)
text(-1.6, 6, "Luzon", pos = 4, cex = 1.1)
text(-1.6, 5, "Visayas", pos = 4, cex = 1.1)
text(-1.6, 4, "Mindanao", pos = 4, cex = 1.1)
text(-1.6, 3, "Difference b/w L & V", pos = 4, cex = 1.1)
text(-1.6, 2, "Difference b/w L & M", pos = 4, cex = 1.1)
text(-1.6, 1, "Difference b/w V & M", pos = 4, cex = 1.1)
axis(1, seq(0, 1, 0.25), labels = NA, lwd = 0.5)
mtext(sprintf("%4.2f", seq(0, 1, 0.25)), side = 1, line = 1, 
      at = seq(0, 1, 0.25), cex = 0.7)
axis(1, seq(1.4, 2.9, 0.25), labels = NA, lwd = 0.5)
mtext(sprintf("%4.2f", seq(-0.5, 1, 0.25)), side = 1, line = 1, 
      at = seq(1.4, 2.9, 0.25), cex = 0.7)
axis(1, seq(3.5, 4.5, 0.25), labels = NA, lwd = 0.5)
mtext(sprintf("%4.2f", seq(0, 1, 0.25)), side = 1, line = 1, 
      at = seq(3.5, 4.5, 0.25), cex = 0.7)
axis(1, seq(4.9, 6.4, 0.25), labels = NA, lwd = 0.5)
mtext(sprintf("%4.2f", seq(-0.5, 1, 0.25)), side = 1, line = 1, 
      at = seq(4.9, 6.4, 0.25), cex = 0.7)
mtext("Percent", side = 1, line = 2.5, at = 0.5, cex = 0.9)
mtext("Percentage points", side = 1, line = 2.5, at = 2.15, cex = 0.9)
mtext("Percent", side = 1, line = 2.5, at = 4, cex = 0.9)
mtext("Percentage points", side = 1, line = 2.5, at = 5.65, cex = 0.9)
mtext("Region", line = 2.5, at = 3.2, font = 2, cex = 1.4)
mtext("Face-to-face", line = 1.5, at = 1.45, font = 2, cex = 1.2)
mtext("Online", line = 1.5, at = 4.95, font = 2, cex = 1.2)
mtext("Approval rate", line = 0.5, at = 0.5, font = 2, cex = 1)
mtext("SDB", line = 0.5, at = 2.15, font = 2, cex = 1)
mtext("Approval rate", line = 0.5, at = 4, font = 2, cex = 1)
mtext("SDB", line = 0.5, at = 5.65, font = 2, cex = 1)
dev.off()

#### cross-tabulation of various subgroup combination (Tables A.2-A.7) ####
## region and socio-economic class
addmargins(table(FtF.data$area, FtF.data$lower), c(1, 2))
round(prop.table(addmargins(table(FtF.data$area, FtF.data$lower), 1), 
                 margin = 1), 3)

## region and internet use
addmargins(table(FtF.data$area, FtF.data$internet.re), c(1, 2))
round(prop.table(addmargins(table(FtF.data$area, FtF.data$internet.re), 1), 
                 margin = 1), 3)

## socio-economic class and Internet use
addmargins(table(FtF.data$lower, FtF.data$internet.re), c(1, 2))
round(prop.table(addmargins(table(FtF.data$lower, FtF.data$internet.re), 1), 
                 margin = 1), 3)

## socio-economic class and perceived neighborhood satisfaction with Duterte
addmargins(table(FtF.data$lower, FtF.data$neighborhood.re), c(1, 2))
round(prop.table(addmargins(table(FtF.data$lower, FtF.data$neighborhood.re), 1), 
                 margin = 1), 3)

## region and perceived neighborhood satisfaction with Duterte
addmargins(table(FtF.data$area, FtF.data$neighborhood.re), c(1, 2))
round(prop.table(addmargins(table(FtF.data$area, FtF.data$neighborhood.re), 1), 
                 margin = 1), 3)

addmargins(table(online.data$area, online.data$neighborhood.re), c(1, 2))
round(prop.table(addmargins(table(online.data$area, online.data$neighborhood.re), 1), 
                 margin = 1), 3)

#### subgroup analyses by ethnicity ####
# significance test of the difference in SDB
ethnicity.result.FtF <- SDB.group.test("list", "nonTag", FtF.data, 
                                       "treatment", "block", "DQ")
round(ethnicity.result.FtF$summary, 3)

FtF.Tagalog.data <- subset(FtF.data, nonTag == 0)
FtF.nonTag.data <- subset(FtF.data, nonTag == 1)

# DQ-based estimate (FtF, Tagalog)
DQ.est.FtF.Tagalog <- mean(FtF.Tagalog.data$DQ, na.rm = TRUE)
DQ.se.FtF.Tagalog <- sqrt(DQ.est.FtF.Tagalog * (1 - DQ.est.FtF.Tagalog) / 
                            sum(! is.na(FtF.Tagalog.data$DQ)))
round(DQ.est.FtF.Tagalog, 3)
round(DQ.est.FtF.Tagalog + qnorm(0.025) * DQ.se.FtF.Tagalog, 3)
round(DQ.est.FtF.Tagalog + qnorm(0.975) * DQ.se.FtF.Tagalog, 3)

# list-based estimate (FtF, Tagalog)
list.result.FtF.Tagalog <- lm_robust(list ~ treatment, data = FtF.Tagalog.data, 
                                     fixed_effects = ~ block)
list.est.FtF.Tagalog <- as.numeric(list.result.FtF.Tagalog$coefficients)
list.se.FtF.Tagalog <- as.numeric(list.result.FtF.Tagalog$std.error)
round(list.est.FtF.Tagalog, 3)
round(list.est.FtF.Tagalog + qnorm(0.025) * list.se.FtF.Tagalog, 3)
round(list.est.FtF.Tagalog + qnorm(0.975) * list.se.FtF.Tagalog, 3)

# SDB estimate (FtF, Tagalog)
SDB.est.FtF.Tagalog <- DQ.est.FtF.Tagalog - list.est.FtF.Tagalog
SDB.se.FtF.Tagalog <- sqrt(list.se.FtF.Tagalog ^ 2 + DQ.se.FtF.Tagalog ^ 2)
round(SDB.est.FtF.Tagalog, 3)
round(SDB.est.FtF.Tagalog + qnorm(0.025) * SDB.se.FtF.Tagalog, 3)
round(SDB.est.FtF.Tagalog + qnorm(0.975) * SDB.se.FtF.Tagalog, 3)
round((1 - pnorm(abs(SDB.est.FtF.Tagalog / SDB.se.FtF.Tagalog))) * 2, 3)

# DQ-based estimate (FtF, nonTag)
DQ.est.FtF.nonTag <- mean(FtF.nonTag.data$DQ, na.rm = TRUE)
DQ.se.FtF.nonTag <- sqrt(DQ.est.FtF.nonTag * (1 - DQ.est.FtF.nonTag) / 
                           sum(! is.na(FtF.nonTag.data$DQ)))
round(DQ.est.FtF.nonTag, 3)
round(DQ.est.FtF.nonTag + qnorm(0.025) * DQ.se.FtF.nonTag, 3)
round(DQ.est.FtF.nonTag + qnorm(0.975) * DQ.se.FtF.nonTag, 3)

# list-based estimate (FtF, nonTag)
list.result.FtF.nonTag <- lm_robust(list ~ treatment, data = FtF.nonTag.data, 
                                    fixed_effects = ~ block)
list.est.FtF.nonTag <- as.numeric(list.result.FtF.nonTag$coefficients)
list.se.FtF.nonTag <- as.numeric(list.result.FtF.nonTag$std.error)
round(list.est.FtF.nonTag, 3)
round(list.est.FtF.nonTag + qnorm(0.025) * list.se.FtF.nonTag, 3)
round(list.est.FtF.nonTag + qnorm(0.975) * list.se.FtF.nonTag, 3)

# SDB estimate (FtF, nonTag)
SDB.est.FtF.nonTag <- DQ.est.FtF.nonTag - list.est.FtF.nonTag
SDB.se.FtF.nonTag <- sqrt(list.se.FtF.nonTag ^ 2 + DQ.se.FtF.nonTag ^ 2)
round(SDB.est.FtF.nonTag, 3)
round(SDB.est.FtF.nonTag + qnorm(0.025) * SDB.se.FtF.nonTag, 3)
round(SDB.est.FtF.nonTag + qnorm(0.975) * SDB.se.FtF.nonTag, 3)
round((1 - pnorm(abs(SDB.est.FtF.nonTag / SDB.se.FtF.nonTag))) * 2, 3)

#### Figure A.11 ####
cairo_pdf("Figure_A11.pdf", width = 4.6, height = 2.6, pointsize = 8)
par(mar = c(4, 0, 3, 1), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-0.7, 1.9), ylim = c(0.5, 3.5), 
     main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = c(seq(0.2, 1, 0.2), seq(1.3, 1.9, 0.2)), col = "gray")
segments(DQ.est.FtF.Tagalog + qnorm(0.025) * DQ.se.FtF.Tagalog, 3.1, 
         DQ.est.FtF.Tagalog + qnorm(0.975) * DQ.se.FtF.Tagalog, 3.1)
points(DQ.est.FtF.Tagalog, 3.1, pch = 15)
text(DQ.est.FtF.Tagalog, 3.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.FtF.Tagalog + qnorm(0.025) * list.se.FtF.Tagalog, 2.9, 
         list.est.FtF.Tagalog + qnorm(0.975) * list.se.FtF.Tagalog, 2.9)
points(list.est.FtF.Tagalog, 2.9, pch = 19)
text(list.est.FtF.Tagalog, 2.9, "List", pos = 1, cex = 0.8)
segments(DQ.est.FtF.nonTag + qnorm(0.025) * DQ.se.FtF.nonTag, 2.1, 
         DQ.est.FtF.nonTag + qnorm(0.975) * DQ.se.FtF.nonTag, 2.1)
points(DQ.est.FtF.nonTag, 2.1, pch = 15)
text(DQ.est.FtF.nonTag, 2.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.FtF.nonTag + qnorm(0.025) * list.se.FtF.nonTag, 1.9, 
         list.est.FtF.nonTag + qnorm(0.975) * list.se.FtF.nonTag, 1.9)
points(list.est.FtF.nonTag, 1.9, pch = 19)
text(list.est.FtF.nonTag, 1.9, "List", pos = 1, cex = 0.8)
segments(SDB.est.FtF.Tagalog + qnorm(0.025) * SDB.se.FtF.Tagalog + 1.3, 3, 
         SDB.est.FtF.Tagalog + qnorm(0.975) * SDB.se.FtF.Tagalog + 1.3, 3)
points(SDB.est.FtF.Tagalog + 1.3, 3, pch = 4)
text(SDB.est.FtF.Tagalog + 1.3, 3, sprintf("%5.3f", SDB.est.FtF.Tagalog), 
     cex = 0.8, pos = 3)
segments(SDB.est.FtF.nonTag + qnorm(0.025) * SDB.se.FtF.nonTag + 1.3, 2, 
         SDB.est.FtF.nonTag + qnorm(0.975) * SDB.se.FtF.nonTag + 1.3, 2)
points(SDB.est.FtF.nonTag + 1.3, 2, pch = 4)
text(SDB.est.FtF.nonTag + 1.3, 2, sprintf("%5.3f", SDB.est.FtF.nonTag), 
     cex = 0.8, pos = 3)
segments(ethnicity.result.FtF$details$SDB.diff.est + qnorm(0.025) * 
           ethnicity.result.FtF$details$SDB.diff.se + 1.3, 1, 
         ethnicity.result.FtF$details$SDB.diff.est + qnorm(0.975) * 
           ethnicity.result.FtF$details$SDB.diff.se + 1.3, 1)
points(ethnicity.result.FtF$details$SDB.diff.est + 1.3, 1, pch = 4)
text(ethnicity.result.FtF$details$SDB.diff.est + 1.3, 1, 
     sprintf("%5.3f", ethnicity.result.FtF$details$SDB.diff.est), 
     cex = 0.8, pos = 3)
text(-0.7, 3, "Tagalog", pos = 4)
text(-0.7, 2, "Non-Tagalog", pos = 4)
text(-0.7, 1, "Difference b/w groups", pos = 4)
axis(1, seq(0.2, 1, 0.2), lwd = 0.5)
axis(1, seq(1.3, 1.9, 0.2), labels = c("0.0", "0.2", "0.4", "0.6"), lwd = 0.5)
mtext("Percent", side = 1, line = 2.5, at = 0.6)
mtext("Percentage points", side = 1, line = 2.5, at = 1.6)
mtext("Approval rate", line = 0.5, at = 0.6, font = 2, cex = 1.1)
mtext("SDB", line = 0.5, at = 1.6, font = 2, cex = 1.1)
dev.off()